1 Introduction

Ultra-relativistic heavy-ion collisions (uRHICs) experiments conducted at the relativistic heavy-ion collider (RHIC) and at large hadron collider (LHC) are the only tools to access experimentally the properties of hot strongly interacting matter. In the last decades it has been reached a general consensus confirmed by several experimental signatures and supported by several theoretical calculations that the matter created in these collisions consists of a strongly interacting quark-gluon plasma (QGP) rather than hadronic matter. One of the most important measurement was the large elliptic flow \(v_2(p_T)= \langle (p_x^2-p_y^2)/(p_x^2+p_y^2) \rangle \) observed [1, 2]. The elliptic flow \(v_2\) is a measure of the anisotropy in momentum space of the emitted particles. It is an observable that encodes information about the Equation of State (EoS) and transport properties of the matter created in these collisions [3,4,5]. The comparison between theoretical calculations and experimental results first within the viscous hydrodynamics framework [6,7,8] and in the recent years also within kinetic transport approach [9,10,11,12,13] have shown that the large value of \(v_2\) is consistent with a matter with a very low shear viscosity to entropy density ratio \(4\pi \eta /s \sim 1-2\). This value is close but larger then the conjectured lower bound for a strongly interacting system, \(\eta /s=1/4\pi \) [14].

In the recent years it has been possible to measure the event-by-event angular distribution of emitted particles. These measurements have made possible to extend this analysis to higher order harmonics \(v_n (p_T)\) [15,16,17]. The origin of these high order harmonics flows is attributed to the fluctuations in the initial geometry [18,19,20,21,22,23,24]. In the recent years most of the theoretical efforts has been focused on the development of the event-by-event viscous hydrodynamics simulations. The comparison between event-by-event viscous hydrodynamical calculations and the experimental results for \(v_n(p_T)\) seems to confirm a finite but not too large average value of \(\eta /s\) with \(4 \pi \langle \eta /s \rangle \sim 1-3\) [21, 22]. Early works, most hydrodynamic simulations of heavy ion collisions, have assumed a constant \(\eta /s\) but a small value of \(\langle \eta /s \rangle \) is not an evidence of the creation of a QGP phase formation. On the other hand, it is known that the \(\eta /s\) of the QGP is expected to have a temperature dependence with a minimum close to the cross over region [25,26,27,28,29]. A phenomenological estimation of its temperature dependence could give further information if the matter created in these collisions undergoes a phase transition [30,31,32].

In the recent years, the CMS, ATLAS, and ALICE collaborations have been able to measure the anisotropic flow coefficients \(\langle v_n \rangle \) in ultra-central heavy ion collisions [17, 33, 34]. These experiments have shown that the triangular flow \(n=3\) appears to be comparable or even larger of the elliptic flow \(n=2\). Such a result poses a problem to hydrodynamical approach that predict a significantly larger \(\langle v_2 \rangle \) w.r.t. \(\langle v_3 \rangle \) in ultra-central collisions. This discrepancy is even more intriguing considering that it occurs for ultra-central collisions where hydrodynamics is expected to work at best. In the recent years, several theoretical calculations within hydrodynamic simulations have shown that in ultra-central collisions the correlations that takes place between the initial eccentricities \(\langle \epsilon _n \rangle \) and the corresponding final anisotropic flows \(\langle v_n \rangle \) are the highest than other centralities. On the other hand, recently within a transport approach it has been pointed-out that at LHC energies and for ultra-central collisions a quite large correlation between \(\langle \epsilon _n \rangle \) and \(\langle v_n \rangle \) similar to viscous hydrodynamics is present up to \(n=5\) at variance with RHIC energies and/or other centralities [24]. At the same time a much larger sensitivity to \(\eta /s(T)\) has been spot. Therefore the study of the anisotropic flows in ultra-central collision offer a powerful tool infer information about the initial geometry of the fireball.

In this paper we study the role of the Equation of State and \(\eta /s (T)\) on the build up of \(v_n(p_T)\) up to the order 5 by using a kinetic transport approach with initial state fluctuations. The paper is organized as follows. In Sect. 2, we introduce the transport approach at fixed shear viscosity to entropy density \(\eta /s\). In Sect. 3, the implementation of the initial state fluctuations in the transport approach. In Sect. 4, we study the role of the EoS on the anisotropic flows and the time evolution of the anisotropic flows \(\langle v_n \rangle \). Finally in Sect. 5 we study the effect of the \(\eta /s(T)\) on the differential \(v_n(p_T)\). In this paper we will show results on \(v_n(p_T)\) for \(n=2,\ldots ,5\) for two different systems: \(Au+Au\) collisions at \(\sqrt{s}=200\) GeV and \(Pb+Pb\) collisions at \(\sqrt{s}=2.76\) TeV at different centralities.

2 Transport approach at fixed \(\eta /s\)

In this work the evolution of the matter created in these collisions is studied employing the kinetic transport theory. This study was performed using a relativistic transport code developed to perform studies of the dynamics of heavy-ion collisions at both RHIC and LHC energies [9, 11, 12, 35,36,37,38,39,40]. Recently this transport code has been extended to include the initial state fluctuations in order to study the elliptic flow \(v_2\) and high order harmonics \(v_n(p_T)\) with \(n > 2\). For a more detailed discussion about the implementation of the initial state fluctuation see [24]. The development of the partonic transport equation has been done only for massless partons [9, 11, 12]. This implies that the matter simulated has an equation of state (EoS) \(\epsilon -3p=0\) which is different to the one evaluated by lQCD [41, 42] which clearly exhibit a large trace anomaly (\(T_{\mu }^{\mu } = \epsilon -3p \ne 0\)) and consequently a sound velocity \(c_{S}^2(T)\) that is significantly smaller than 1 / 3 also in the initial stages of a Heavy Ion Collisions reaching about 1 / 10 at \(T \approx T_C=155 \, MeV\). In this paper we extend the solution of the Relativistic Boltzmann Transport (RBT) equation to massive partons which allows to simulate a fluid with an EoS similar to the recent lQCD calculations.

Therefore, the evolution of the phase-space distribution function f(xpt) is given by solving the following RBT equation:

$$\begin{aligned}&p^{\mu }\, \partial _{\mu } f(x,p)= C[f] + S[f_0] \end{aligned}$$
(1)

with \(p^{\mu }p_{\mu }=m^2\) and with \(f_0\) the initial distribution function (that can be of non-equilibrium). While C[f] is the Boltzmann-like collision integral. In this paper we have considered only the \(2 \leftrightarrow 2\) processes. The C[f] can be written as,

$$\begin{aligned} C[f]=\int \frac{d^3{p}_{2}}{E_{2}} d\varOmega \frac{1}{2}\sqrt{s(s-4m^2)} \sigma (s,\varTheta ) (f_{1^\prime }f_{2^\prime } - f f_2)\nonumber \\ \end{aligned}$$
(2)

where \(\sigma (s,\varTheta )\) is the differential cross section which is related to the total cross section by \(\sigma _{tot}=\int d\varOmega \sigma (s,\varTheta )\). The transport theory permits to study the effects of the microscopical processes fixed by matrix element or cross sections on the observables. This is however not the aim of this work where instead we determine an effective cross section \(\sigma _{tot}\) as a tool to fix the \(\eta /s\) of the system. This approach is inspired by the success of the hydrodynamical approach that has shown the key role played by the \(\eta /s\) to describe the experimental data [21, 22, 43]. Therefore we use the RBT equation as an approach converging to hydrodynamics for small scattering relaxation time \(\tau \sim 1/\sigma \rho \) (or small \(\eta /s\)) [44,45,46]. On the other hand, it permits to study directly the impact of a T dependence of \(\eta /s\) on observables like the anisotropic flows \(v_2(p_T)\) and \(v_3(p_T)\) which is the main focus of this paper. However, this approach it is not based on an ansatz for the viscous corrections for the phase-space distribution function \(\delta f\) and it is not affected by limitation in the transverse momentum range in order to ensure that \(\delta f/f<< 1\). Also the RBT approach is naturally valid also at large \(\eta /s\) or for \(p_T>> T\) and it is possible to include initial non equilibrium effects to the initial distribution function f(xp) [37, 38, 47]. Furthermore, as will be discussed below, the kinetic freeze-out can be determined self-consistently by increasing \(\eta /s(T)\) at low temperature which gives a smooth switching-off of the scattering rates. The current disadvantage of the present approach is that hadronization has not yet been included and it will be the main topics of forthcoming studies. In the transport calculations shown in this paper we use massive particles providing a soften equation of state respect to the massless case with a decreasing speed of sound when the crossover region is approached similar to the one of lQCD calculations [41, 42].

To study the expansion dynamics of the fireball with a certain \(\eta /s(T)\), we determine locally in space and time the total cross section \(\sigma _{tot}\) according to the Chapmann–Enskog theory [48, 49]. For the results shown in this paper we have considered an isotropic cross section. As shown in [48,49,50] for a massive gas in the Chapmann–Enskog approximation at the first order \(\eta \) is given by

$$\begin{aligned} \eta = f(z) \frac{T}{\sigma _{tot}} \end{aligned}$$
(3)

with \(z=m/T\) while the function f(z) is given by

$$\begin{aligned} f(z) = \frac{15}{16}\frac{\left[ z^2 K_3(z)\right] ^{2}}{(15z^2+2)K_2(2z)+(3z^3+49z)K_3(2z)} \end{aligned}$$
(4)

where \(K_n\)-s are the modified Bessel functions. The entropy density for a massive system is given by \(s=\rho \big ( 4 + z K_1(z)/K_2(z) \big )\). Notice that in the ultra-relativistic limit \(z\rightarrow 0\) the \(f(z) \rightarrow 1.2\) and we recover the massless limit for the \(\eta \). As shown in [36] the expression for \(\eta \) in Eq. (3) is in quite good agreement at level of \(3-5\%\) with the Green-Kubo formula and it describes correctly the \(\eta /s\) of the system in the range of interest of HIC. Notice that in this limit the \(\eta /s\) ratio coincides with the one obtained in the relaxation time approximation.

Fig. 1
figure 1

Different parametrizations for \(\eta /s(T)\). The orange area refers to the quasi-particle model predictions for \(\eta /s\) [51]. The three different lines indicate different possible T dependencies used in this paper. Symbols are as in the legend. See the text for more details

In Fig.1 it is shown a collection of different theoretical calculations about the T dependence of the \(\eta /s\) ratio. In literature there are several indications that \(\eta /s\) should have a particular behaviour with the temperature [25,26,27, 30, 31, 52]. As shown in Fig.1 in general \(\eta /s\) should have a typical behaviour of phase transition with a minimum close to the critical temperature \(T_C\) [30,31,32]. Estimations of \(\eta /s\) within chiral perturbation theory for a meson gas [25, 26], have shown that at temperature lower than the critical temperature the \(\eta /s\) is a decreasing function with the temperature, see down-triangles in Fig.1. Similar results have been obtained from extrapolation of heavy-ion collisions data at intermediate energies, see HIC-IE diamonds in Fig.1. At higher temperature \(T>T_c\) lQCD calculation on pure gauge have shown that in general \(\eta /s\) becomes an increasing function with the temperature [27,28,29], see up-triangles, circles and squares in Fig.1. However due to the large error bars in the lQCD results for \(\eta /s\) at moment it is not possible to infer a clear temperature dependence in the QGP phase.

In Fig.1 we show the three different T dependencies of \(\eta /s\) studied in this paper. In particular we have considered the following cases: one of constant \(4\pi \eta /s=1\) during all the evolution of the system shown by dot dashed line in Fig.1 another one of \(4\pi \eta /s=1\) at higher temperature in the QGP phase and an increasing \(\eta /s\) in the cross over region towards the estimated value for hadronic matter \(4\pi \eta /s \approx 6\) [26, 53] and shown by solid line. The third one is shown in Fig.1 by the dot dashed line. In this case we consider the increase of \(\eta /s\) at higher temperature with a linear temperature dependence and a minimum close to the critical temperature with a temperature dependence similar to that expected from general considerations as suggested by quasi particle models [51, 54] or by recent calculations on pure gauge [29]. Notice that in our approach an increase of \(\eta /s\) at lower temperature \(0.8 T_C \le T \le 1.2 T_C\) permits to account for a smooth realistic kinetic freeze-out (f.o.) because at lower T the mean free path \(\lambda \) goes like \(\lambda \propto \frac{\eta }{s} \frac{1}{T}\) and the total cross section goes like \(\sigma \propto (\eta /s)^{-1}\). In the following discussion the term f.o. means that we take into account the increase of \(\eta /s\) at \(T < 1.2 T_C\).

3 Initial conditions

In the following paper we will consider two systems at different centralities: \(Au+Au\) collisions at \(\sqrt{s_{NN}}=200\) GeV produced at RHIC and \(Pb+Pb\) collisions at \(\sqrt{s_{NN}}=2.76\) TeV at LHC. In this section we discuss the implementation of the initial state fluctuations in the transport approach, for more details about the implementation of the initial fluctuations in our approach see Ref. [24]. To generate an event by event initial profile we use the Monte-Carlo Glauber model. For a fixed centrality class we determine the average impact parameter by the Monte Carlo Glauber model. We use this to generate event by event \(N_{part}\) according to the Monte Carlo Glauber model. We employ the geometrical method to determine the the initial number of participant nucleons. Within this method two nucleons collide each other if the relative distance in the transverse plane is \(d_{T} \le \sqrt{\sigma _{NN}/\pi }\) where \(\sigma _{NN}\) is the nucleon-nucleon cross section. In our calculation \(\sigma _{NN}=4.2 \, fm^2\) is taken for RHIC at \(\sqrt{S_{NN}}=200\) GeV and \(\sigma _{NN}=7.0 \, fm^2\) for LHC at \(\sqrt{s_{NN}}=2.76\) TeV. The resulting initial nucleon distribution is converted in a smooth one by assuming a Gaussian distribution centered in the nucleon position as done in several hydrodynamics calculations [8, 19, 20, 43, 55]. In our calculations we convert the information of the nucleon distribution into parton density in the transverse plane \(\rho _T(x,y)\) by the following expression

$$\begin{aligned} \rho _T(x,y) = k \sum _{i=1}^{N_{part}} \exp {\bigg [-\frac{(x-x_i)^2+(y-y_i)^2}{2\sigma _{xy}^2}\bigg ]} \end{aligned}$$
(5)

where \(N_{part}\) is the number of participants and k is an overall normalization factor fixed in order to reproduce the final hadron multiplicity \(dN_{ch}/dy\). In the results shown in this paper we have fixed the Gaussian width of the fluctuations to \(\sigma _{xy} = 0.5 \, fm\). In our calculation we assume a longitudinal boost invariant distribution from \(y=-2.5\) to \(y=2.5\).

Fig. 2
figure 2

Initial eccentricities \(\langle \epsilon _n \rangle \) for \(n=2,\ldots ,5\) as a function of the impact parameter. Different symbols are for different harmonics n. In the inset panel it is shown the \(\langle \epsilon _n \rangle \) as a function of n for ultra-central collisions

Fig. 3
figure 3

Initial temperature in the transverse plane T(xy) at mid rapidity for two typical events for \(Au+Au\) at \(\sqrt{s_{NN}}=200\) GeV (upper panel) and \(Pb+Pb\) at \(\sqrt{s_{NN}}=2.76\) TeV (lower panel). These plots are for central collision

The transverse density \(\rho _T(x,y)\) fixes in each event the initial anisotropy in coordinate space that is quantified in terms of the following coefficients \(\epsilon _n\) [19, 20]:

$$\begin{aligned} \epsilon _{n}=\frac{\sqrt{\langle r_T^n \cos {(n \phi )} \rangle ^2 + \langle r_T^n \sin {(n \phi )} \rangle ^2}}{\langle r_T^n\rangle } \end{aligned}$$
(6)

where \(r_T=\sqrt{x^2+y^2}\) and \(\phi =arctan(y/x)\) is the polar coordinate in the transverse plane. In Fig. 2 it is shown the initial spatial anisotropies \(\langle \epsilon _n \rangle \) for \(n=2,\ldots ,5\) as a function of the impact parameter. In the inset panel it is shown \(\langle \epsilon _n \rangle \) as a function of the order of the harmonic for central collision. The second coefficient \(\langle \epsilon _2 \rangle \) shows a stronger dependence with the impact parameter with respect to the other coefficients and it dominates for large impact parameter. This is due to the global almond shape of the fireball while the other harmonics have most of their origin in the fluctuations of the positions of the nucleons. For more central collisions \(b \le 2.5\) fm we observe an opposite behaviour and the \(\langle \epsilon _2 \rangle \) is even smaller than the other harmonics. This is due to the fact that the effect of the elliptic overlap region disappears it becomes more difficult to have fluctuations of the positions of the nucleons along only one preferential direction. In particular, as shown in the inset panel, for central collision larger is the order of the harmonics larger is the corresponding initial eccentricities. However this is a general behaviour for the eccentricities for central collisions and similar ordering happens also for different bin selections. For the initialization in momentum space at RHIC (LHC) energies we have considered for partons with transverse momentum \(p_T \le p_0=2\) GeV (3 GeV) a thermalized spectrum in the transverse plane. Assuming the local equilibrium the initial local temperature in the transverse plane T(xy) is evaluated by using the standard thermodynamical relation for massive particles \(\rho _T(x,y)=\gamma T^3 z^2 K_{2}(z)/(2 \pi ^{2})\) with \(z=m/T\) and \(\gamma =2 \times (N_c^2-1) + 2 \times 2 \times N_c \times N_f=40\) with \(N_c=3\) and \(N_f=2\). In Fig. 3 it is shown the initial temperature profile in the transverse plane for a given event and for \(Au+Au\) collisions at \(\sqrt{s_{NN}}=200\) GeV (upper panel) at initial time of \(\tau _0=0.6\) fm/c and for \(Pb+Pb\) collisions at \(\sqrt{s_{NN}}=2.76\) TeV (lower panel) at initial time \(\tau _0=0.3\) fm/c. These plots are at mid rapidity and for central collisions.

While for partons with \(p_T > p_0\) we have assumed the spectrum of non-quenched minijets according to standard NLO-pQCD calculations with a power law shape [56, 57]. The initial transverse momentum of the particles is distributed uniformly in the azimuthal angle. We fix the starting time of the simulation to \(\tau _0=0.6\) fm/c for RHIC and \(\tau _0=0.3\) fm/c for LHC as commonly done also in hydrodynamical approaches.

In our simulations we have used \(N_{event}=1000\) events for each centrality class. This number is enough to get stable results for the spectra, differential elliptic flow and high order flow coefficients \(v_n(p_t)\) for \(n=2,\ldots ,5\). Moreover, in the calculation shown in this paper have checked the convergence of our results for \(v_n\) with the lattice spacing of the calculation grid and the number of the test particles \(N_{test}\). In the results hown in this paper we have used a calculation grid with a transverse area of the cell \(A_T=0.12\, \mathrm{{fm}}^2\) and the total number of test particles is fixed to be \(N_{test}=2 \cdot 10^6\) for each event.

Fig. 4
figure 4

\(p/\epsilon \) ratio as a function of the energy density. The orange band are the lQCD data taken from the Wuppertal–Budapest collaboration [58]. The solid green and dashed blue lines are respectively for LHC and RHIC energies for central collisions and mid rapidity

4 Role of the equation of state on the \(v_n\)

As well known the anisotropic flow coefficients \(v_{n}(p_T)\) are observable sensitive to the bulk properties of the created matter and they encode information about the EoS and transport properties of the medium like for example gives information about \(\eta /s\) ratio. In this section we discuss the role of the EoS on the build up of the \(v_{n}(p_T)\). To extract information about the EoS we evaluate the average energy-momentum tensor.

$$\begin{aligned} \langle T^{\mu \nu } \rangle = \frac{1}{V N_{test}}\sum _{i=1}^{N_{tot}} \frac{p_{i}^{\mu } p_{i}^{\nu }}{p_{i}^0} \end{aligned}$$
(7)

We evaluate the energy-momentum tensor in a cylinder of radius \(r_T=1.5\) fm and in longitudinal direction \(|\eta |<0.3\). We evaluate the energy density by \(\epsilon = \langle T^{00} \rangle \) and the pressure by assuming that \(p= \frac{1}{3} \sum _{i} \langle T^{ii} \rangle \). In Fig. 4 it is shown the comparison of \(p/\epsilon \) ratio obtained dynamically in our simulations with the one of lQCD for both RHIC energies in \(Au+Au\) collisions at \(\sqrt{s_{NN}}=200\) GeV and at LHC energies in \(Pb+Pb\) collisions at \(\sqrt{s_{NN}}=2.76\) TeV for central collisions. The results shown a reasonably agreement with lQCD data for both energies. In these calculations we have fixed \(m=0.5\) GeV that are in agreement with the masses extracted within quasi-particle models [51, 59]. As shown, this implies that in our simulations dynamically we describe a systems that approximatively have the equation of state as the one in lQCD calculations. The initial fluctuations in density induce fluctuations in energy density and the corresponding value of \(p/\epsilon \) can fluctuate event-by-event. In Fig. 4 the blue and green bands show the possible range of trajectories covered in our simulations. The band is due to the fact that locally the relation between pressure and energy density can be modified by non equilibrium effects.

In this paper the anisotropic flows \(v_n(p_T)\) with \(n=2,\ldots ,5\) have been calculated using the event plane method as

$$\begin{aligned} v_n=\langle \cos {[n(\phi -\varPsi _n)]}\rangle \end{aligned}$$
(8)

where the momentum space angles \(\varPsi _n\) are given by

$$\begin{aligned} \varPsi _n=\frac{1}{n}\arctan {\frac{\langle \sin {(n\phi )}\rangle }{\langle \cos {(n\phi )}\rangle }} \end{aligned}$$
(9)

while the eccentricities are evaluated using Eq. (6) as done in [8, 19, 20].

As a first step we study the time evolution of the initial eccentricities \(\langle \epsilon _n \rangle \) and the corresponding integrated anisotropic flows \(\langle v_n \rangle \). Where the average \(\langle ... \rangle \) is over all the events.

Fig. 5
figure 5

Results for \(Pb+Pb\) collisions at \(\sqrt{s_{NN}}=2.76\) TeV for (20–30)% centrality collisions and for the case \(4\pi \eta /s=1\) and the kinetic freeze out (f.o.). Left panel: time evolution of initial eccentricities \(\langle \epsilon _n \rangle \) different symbols are for different harmonics solid lines are for the massive case while dashed lines are for massless case. Right panel: time evolution of \(\langle v_n\rangle \) at mid rapidity respectively for (20–30)% centrality collisions. Same legend as in left panel

In the left panel of Fig. 5 we plot the time evolution of the initial anisotropies in coordinate space \(\langle \epsilon _{n} \rangle \) for both cases massless (dashed lines) and massive case (solid lines). In the right panel of Fig. 5 instead we show the time evolution for the corresponding final anisotropic flow coefficients \(\langle v_{n}\rangle \). These calculation are for LHC energies at mid rapidity and for (20–30)% centrality collisions. Results show that in general for both massless and massive case the \(v_n(p_T)\) are built-up within the first few fm/c and for both cases different harmonics have a different formation time with the following ordering \(t_{v_2}<t_{v_3}<t_{v_4}\cdots \) as shown in Ref. [24], but here we consider the impact of the EoS. Moreover as shown in the right panel of Fig. 5 we observe that for massless case the system is more efficient to convert the initial anisotropy in coordinate space into final anisotropy in momentum space where the effect of a realistic EoS is to reduce the \(\langle v_{n}\rangle \). This effect increase with the increasing of the order of the harmonics. The larger effect for larger n is related to the fact that higher harmonics have larger formation time and therefore they develop in a region where the \(p/\epsilon \) or the speed of sound is smaller. In fact for example the second coefficient \(\langle v_{2}\rangle \) is almost completely developed within 3–4 fm/c where the \(p/\epsilon \approx 0.3\) and close to the conformal limit and the effect increases for large n because higher harmonics develop later at lower temperature where \(p/\epsilon < 0.3\).

Fig. 6
figure 6

Differential \(v_n(p_T)\) for \(Pb+Pb\) collisions at \(\sqrt{s_{NN}}=2.76\) TeV for (20–30)% centrality collisions and for the case \(4\pi \eta /s=1\) and the kinetic freeze out (f.o.). Different symbols are for different harmonics n. Solid lines are for the massive case while the dashed lines are for the massless case

In Fig. 6 we compare the final \(v_n(p_T)\) for both massless (dashed lines) and massive case (solid lines). These results are for \(Pb+Pb\) collisions at \(\sqrt{s_{NN}}=2.76\) TeV and for (20–30)% centrality collisions. We observe that for the elliptic flow \(v_2\) a mass ordering at low transverse momentum where larger is the mass smaller is the corresponding elliptic flow which is typical of hydrodynamic expansion that boost strongly massive particles [4, 60]. The explanation of the mass ordering is due to the fact that at low \(p_T\) the elliptic flow is \(v_{2}(p_T) \propto p_T -\langle \beta _{T}\rangle m_T\) therefore for light particles \(m_T = \sqrt{p_T^2+m^2} \approx p_T\) and \(v_2(p_T) \propto p_T\) while for massive particle the coupling of \(m_T\) with the radial flow gives at the same \(p_{T}\) a smaller \(v_2(p_T)\) and it is translated in a non linear behaviour at low \(p_T\). Moreover a similar mass ordering it is shown also for higher harmonics \(v_n(p_T)\). At low \(p_T\) the effect is to give a non linear behaviour for the second and third coefficients. Moreover we observe that the sensitivity to the EoS of the \(v_n(p_T)\) increases with the increasing of the harmonic n where for example the reduction is less of 10% for the third harmonic and about 15% for \(v_{4}(p_T)\) and 25 % for \(v_{5}(p_T)\). The enhancement of the sensitivity with the increase of the order of the harmonic is again related to the different formation time where for large n the corresponding \(v_n(p_T)\) develops later where the speed of sound \(c_{s}^2\) is smaller compared to the one in the conformal limit with \(p/\epsilon \approx 0.25\) at \(t \approx \, 3 fm/c\), see Fig. 4.

5 Effects of \(\eta /s(T)\) on the \(v_n(p_T)\)

Differential anisotropic flow coefficients \(v_n(p_T)\) are observables that carry out more information about the fireball created in the heavy ion collisions. Moreover they are observables sensitive to the transport properties of the medium like the \(\eta /s\) ratio. In the following discussion we will study the effect of the \(\eta /s\) on the build up of the elliptic flow \(v_2(p_T)\) and on the high order harmonics \(v_3(p_T)\), \(v_4(p_T)\) and \(v_5(p_T)\). In particular we study two different systems the one created at RHIC for \(Au+Au\) collisions at \(\sqrt{s} = 200\) GeV and the one at LHC energies for \(Pb+Pb\) at \(\sqrt{s}= 2.76\) TeV. The following calculations are for the massive case which gives a reasonable agreement with the lQCD EoS, see Fig. 4.

Fig. 7
figure 7

Differential \(v_n(p_T)\) at mid rapidity and for (20–30)% collision centrality. The comparison is between the two systems: \(Au+Au\) at \(\sqrt{s}=200\) GeV (left panel) and \(Pb+Pb\) at \(\sqrt{s}=2.76\) TeV (right panel). Different colors are for different harmonics n. The dot dashed lines refer to the case with a constant \(\eta /s=(4 \pi )^{-1}\) during all the evolution. The solid lines refer to the case with \(\eta /s=(4 \pi )^{-1}\) at higher temperature and kinetic f.o. (see Fig. 1) while the dashed lines refer to the case with \(\eta /s \propto T\) at higher temperature and the kinetic f.o. at lower temperature

In Fig. 7 we show the effect of the \(\eta /s(T)\) on the \(v_{n}(p_T)\) at mid rapidity and for (20–30)% centrality for both RHIC energies for \(Au+Au\) collisions at \(\sqrt{s} = 200\) GeV (left panel) and LHC energies for \(Pb+Pb\) collisions at \(\sqrt{s}= 2.76\) TeV (right panel). Similarly to what has been obtained in viscous hydrodynamical calculations or in transport simulation the increase of the viscosity of the medium has the effect to reduce the produced \(v_{n}(p_T)\). These results are qualitatively in agreement with previous results within kinetic transport approach in the conformal limit [24]. Comparing the dot dashed lines with the solid lines in the left panel of Fig. 7, we observe that at RHIC energies the \(v_n(p_T)\) are sensitive to the increase of the \(\eta /s\) at lower temperature close to the cross over region. As one should expect by the different formation time, this sensitivity increase with the order of the harmonics. On the other hand at LHC energies, right panel of Fig. 7, the \(v_n(p_T)\) are more sensitive to the increase of \(\eta /s\) at high temperature as shown by a larger suppression of the \(v_n(p_T)\) for the case where \(\eta /s \propto T\) in the QGP phase (dashed lines).

Fig. 8
figure 8

Differential \(v_n(p_T)\) at mid rapidity and for ultra-central collisions. The comparison is between the two systems: \(Au+Au\) at \(\sqrt{s}=200\) GeV (left panel) and \(Pb+Pb\) at \(\sqrt{s}=2.76\) TeV (right panel). Different colors are for different harmonics n. The dot dashed lines refer to the case with a constant \(\eta /s=(4 \pi )^{-1}\) during all the evolution. The solid lines refer to the case with \(\eta /s=(4 \pi )^{-1}\) at higher temperature and kinetic f.o. while the dashed lines refer to the case with \(\eta /s \propto T\) at higher temperature and the kinetic f.o. at lower temperature

In the recent years, it has been possible to access experimentally to the ultra-central collisions. In ultra-central collisions the initial asymmetry in coordinate space measured by \(\epsilon _n\) comes from the fluctuations in the initial geometry and there is no coupling with the global overlap shape. Moreover, it has been shown within event-by-event ideal and viscous hydrodynamics calculations [55, 61, 62] that for these collision centralities the final integrated anisotropic flows coefficient \(\langle v_{n}\rangle \) are strongly correlated to the initial eccentricities \(\langle \epsilon _{n}\rangle \) with a linear correlation coefficient \(C(\epsilon _n,v_n)\approx 1\). Similar results have been obtained also within event-by-event transport approach [24].

In Fig. 8 it is shown the comparison of \(v_{n}(p_T)\) produced at RHIC energies in \(Au+Au\) collisions at \(\sqrt{s}=200\) GeV (left panel) with the one produced at LHC energies in \(Pb+Pb\) collisions at \(\sqrt{s}=2.76\) TeV (right panel) for ultra-central collisions. Comparing Fig. 7 with Fig. 8 one can observe that at low \(p_T\) both centralities the \(v_{n}(p_T) \propto p_{T}^{n}\) while at hight \(p_T\) for ultra-central collisions we observe that the elliptic flow \(v_{2}(p_T)\) shows a saturation while \(v_{n}(p_T)\) for \(n \ge 3\) have an almost linear \(p_T\) dependence. These results are in qualitative agreement with what has been observed experimentally. Similarly to the results at mid peripheral collisions at RHIC energies the \(v_n(p_T)\) are more sensitive to the value of \(\eta /s\) at lower temperature. Moreover, we observe that for central collisions the \(v_n(p_T)\) are more sensitive with respect to more peripheral collisions with a suppression of about 30–35% for \(n > 3\) in central collision and of about 20–25% for mid peripheral collisions. On contrary at LHC energies we observe that reduction of \(v_{n}(p_T)\) due to the increase of \(\eta /s\) in the QGP phase (dashed lines) is strongly enhanced for ultra-central collisions. In particular for \(n \ge 3\) the reduction for central collisions is about 30–35% against a reduction of about 10% for mid peripheral collisions. It is indeed remarkable notice that the 30% effect is determined by a slowly linear rising of \(\eta /s\) with T and with a minimum close \(T_C\) of \(\eta /s\approx 1/(4\pi )\) as the dashed lines in Fig. 1. Moreover in central collisions at LHC energies higher harmonics acquire a larger sensitivity to the value of the viscosity in the QGP phase. Such a sensitivity is not present at RHIC energy and at LHC in non-central collisions. On the other hand this does not appear for the \(\langle \epsilon _n \rangle \) but only at \(p_T \gtrsim 1\) GeV.

Fig. 9
figure 9

\(\langle v_{n} \rangle \) as a function of the order of the harmonic n for ultra-central collisions and for \(Au+Au\) collisions at \(\sqrt{s_{NN}}=200\) GeV (left panel) and \(Pb+Pb\) collisions at \(\sqrt{s_{NN}}=2.76\) TeV (right panel). Different symbols are for different T dependence of the \(\eta /s\) ratio: black circles are for \(4\pi \eta /s=1\) during all the evolution of the system, red squares for the case \(4\pi \eta /s=1\) and the kinetic f.o. while green triangles for \(\eta /s \propto T\) at high temperature and the kinetic f.o. at low temperature. The orange diamonds are the results for the massless case with \(4\pi \eta /s=1\) and the kinetic f.o. at lower temperature

The comparison between RHIC and LHC highlights the role played by the EoS and the role of the collision energy in the build-up of \(\langle v_{n} \rangle \) in ultra central collisions. In ultra-central collisions, anisotropic flows coefficients are generated by density fluctuations in the initial state rather than by geometric overlap effects. On the other hand, the predicted anisotropic flow coefficients for charged hadrons with viscous hydrodynamics calculations with both the Monte Carlo Glauber (MC-Glb) and Monte Carlo Kharzeev–Levin–Nardi (MC-KLN) models produce initial fluctuations that gives final anisotropic flow spectrum as a decreasing function with the order of the harmonic n. On the other hand the CMS Collaboration has measured the anisotropic flow coefficients of charged hadrons in ultra-central Pb+Pb collisions at \(\sqrt{s}=2.76\) TeV corresponding to 0–0.2% centrality. The experimental measurement have shown that the triangular flow \(n=3\) appears to be comparable or even larger of the elliptic flow \(n=2\). Such a result poses a problem to hydrodynamical approach that predict a larger \(\langle v_2 \rangle \) w.r.t. \(\langle v_3 \rangle \) in ultra-central collisions. In Fig. 9 it is shown \(\langle v_{n} \rangle \) as a function of the order of the harmonic n for ultra-central collisions and for \(Au+Au\) collisions at \(\sqrt{s_{NN}}=200\) GeV (left panel) and \(Pb+Pb\) collisions at \(\sqrt{s_{NN}}=2.76\) TeV (right panel) As shown in the right panel of Fig. 9 by comparing circles with diamonds the effect of a soften Equation of State is to reduce the build-up of the \(\langle v_{n} \rangle \). In particular we observe a greater suppression for higher harmonics where for example for the massless case we get that \(\langle v_{4} \rangle \approx \langle v_{2} \rangle \) while for the massive case \(\langle v_{4} \rangle < \langle v_{2} \rangle \) with a suppression of about 15%. On the other hand the role of the collision energy is to reduce the production of the anisotropic flows and moreover it can play a role on the peak for \(n=3\). The different behaviour between RHIC and LHC is due to the fact that the degree of correlation between the initial \(\langle \epsilon _n \rangle \) and the final anisotropic flows produced \(\langle v_{n} \rangle \) is different. In fact the linear correlation coefficient between them \(C(\epsilon _n,v_m)\) is a decreasing function with the collision energy and impact parameter as shown for massless partons in [24]. This implies that for ultra-central collisions and at LHC energies the degree of correlation is maximum with \(C(\epsilon _n,v_m)\approx 1\) for \(n=2,\ldots ,4\). Therefore at higher energies like at LHC energies the \(\langle v_n \rangle \) keep more information about the initial anisotropies in the coordinate space and in particular they tend to keep the initial ordering of the eccentricities where for ultra-central collisions we have \(\langle \epsilon _2 \rangle< \langle \epsilon _3 \rangle < \langle \epsilon _4 \rangle \cdots \), see Fig. 2. At LHC it has been observed a \(\langle v_{2} \rangle \lesssim \langle v_{3} \rangle \) for ultra-central collisions. Our simulation show that such a feature should disappear at RHIC because of the impact of the increasing \(\eta /s(T)\) in the cross-over region. On the other hand at LHC a \(\langle v_{2} \rangle \lesssim \langle v_{3} \rangle \) survives because the freeze-out dynamics implied by \(\eta /s(T)\) does not affect the \(\langle v_{n} \rangle \) at high energy.

6 Conclusions

In this paper we have investigated the build up of the anisotropic flows \(v_{n}(p_T)\) for \(n=2,3,4\) and 5 within an event-by-event transport approach. We have studied the role of the EoS, by using a finite partonic mass, on the time evolution of \(\langle v_n \rangle \) and the initial eccentricities \(\langle \epsilon _n \rangle \) and on the differential \(v_{n}(p_T)\). We have found that the effect of the mass is to give a non linear behaviour at low \(p_T\) for the elliptic and triangular flow. In general the mass has the effect to reduce the final \(v_{n}(p_T)\) produced. The system is less efficient in converting the initial anisotropy in coordinate space into final anisotropic flows. Due to the different formation time of the harmonics the efficiency is lower for higher harmonics that start to develop at smaller speed of sound. Moreover, we have also studied the effect the temperature dependence of the \(\eta /s\) ratio on \(v_{n}(p_T)\) for two different beam energies: at RHIC for \(Au+Au\) collisions at \(\sqrt{s}=200\) GeV and at LHC for \(Pb+Pb\) collisions at \(\sqrt{s}=2.76\) TeV. We have found that at RHIC the \(v_n(p_T)\) are more affected by the value of \(\eta /s\) in the cross over region (\(T<1.2 T_C\)) and the sensitivity increases with the order of the harmonics. At LHC energies we get that almost all the \(v_n(p_T)\) develop in the QGP phase and they have less contamination of the value of \(\eta /s\) in the cross-over region and we observe a weaker sensitivity to the T dependence of the \(\eta /s\). These results are in qualitative agreement with the results obtained for the massless case. Such a scenario changes for ultra-central collisions. In particular in these collision an enhancement of the sensitivity to the value of \(\eta /s\) at lower temperature for RHIC is observed while at LHC the \(v_n(p_T)\) are more sensitive to the value at high temperature. In general we found an enhancement of the sensitivity of the \(v_n(p_T)\) for \(n=2,\ldots ,5\) that can reaches about a 30–35%.