Using the non-hydrodynamic mode to study the onset of hydrodynamic behavior in ultraperipheral symmetric nuclear collisions

With the attempts of extending the hydrodynamic framework of heavy-ion collision to proton-proton and other small and low energy systems, we are confronted with the question of how small the system can get and still be safely modelled as a fluid. One of the transport coefficients required in the $2^{nd}$ order relativistic viscous hydrodynamics is the shear relaxation time, inclusion of which solves the causality violation problem in the Navier-Stokes equation. In phenomenological studies this coefficient has been taken as a constant and much attention has gone into finding and fixing the shear viscosity to entropy density ratio, $\eta/s$. This transport coefficient also happens to control the non-hydrodynamic mode of the out-of-equilibrium hydrodynamics theory. It has been predicted that for decreasing system size, observables become sensitive to variation in shear relaxation time as a result of increasing dominance of non-hydrodynamic mode, which could potentially indicate breakdown of hydrodynamics. In this study, we try to test this prediction in the peripheral Pb-Pb collisions at $2.76$ TeV and Au-Au collisions at $200$ GeV, with IPGlasma initial condition and $(2+1)-$Dimensional viscous hydrodynamics. We find that elliptic flow does show adequate sensitivity to variation in relaxation time for decreasing system size. The multiplicity rapidity density limit for applicability of hydrodynamics is found to be around $dN/dy\approx10$, with the possibility of refinement in this value given a way to improve the centrality resolution in experimental data for referencing in peripheral collisions.


INTRODUCTION
The fact that baryons have internal structure directly leads to the notion that a bulk medium of sub-nucleonic degrees of freedom should exist [1, 2]. An energy density of about 0.7 GeV/fm 3 is required to free up quarks from the nucleons [3,4]. We now have convincing signs from experiments at the BNL Relativistic Heavy Ion Collider (RHIC) and the CERN Large Hadron Collider(LHC) that indicate a deconfined state of quarks and gluons called quark-gluon plasma (QGP) is formed for a sufficient distinguishable duration. Low-order hydrodynamic constitutive relations apparently explain the experimental observables of such a dynamic system quite well, even though there is a sizable pressure anisotropy. This applicability of low-order hydrodynamics has been referred to as hydrodynamization, 1 , to distinguish it from local thermalization [5,6]. Experimental confirmation of strangeness enhancement [7], elliptic flow [8] and jet quenching [9,10] as the early indicators was subsequently followed by confirmation of other signatures like quarkonia suppression. Efforts now are directed towards quantitatively fixing the boundaries of various regions of Quantum Chromodynamics(QCD) phase diagram [11] and deducing the properties of QGP [12].
There are challenges involved in analytically solving non-perturbative QCD making the proof of deconfinement intractable [13]. Hence the progress in modelling a medium of quarks and gluons from first principles has been limited. Lattice QCD, even though computationally intensive, has been of help in understanding deconfinement, and other low density phenomena where the numerical sign problem does not affect the calculations [14,15]. For now, phenomenological models aided by lattice QCD seem to be the right approach in modelling such a complex system. The use of hydrodynamics in modelling the transient QGP stage has been quite surprising [16]. However, hydrodynamics as an effective theory for heavy-ion collisions, has evolved tremendously, especially in the last two decades. For an in-depth review of the hydrodynamics in heavy-ion collisions, please lookup Refs. [17][18][19][20][21][22][23]. Apart from the traditional conversation equation approach, hydrodynamics can also be derived as a microscopic theory in the limit e.g., starting from kinetic theory or any QFT like QCD provided its dynamics show a quasi-universality at a large time scale [21]. This microscopic theory approach also helps in fixing the transport coefficients of the theory [18].
The energy momentum tensor for such a theory in a non-equilibrium state is decomposed as; where the first and second term represents the equilibrium state of T µν and the deviation from the equilibrium, respectively. Under linear response theory, the second term can be expanded as; where ω is the angular frequency and, k is the momentum, has singularities. The solution of the δ T µν (x) integral at late times has a contribution in terms of complex singular frequency in the ω-plane; where, ω h is the real part of frequency at singularity corresponding to excitation of equilibrium plasma, also called hydrodynamic mode frequency. ω nh is termed as transient mode or non-hydrodynamic mode frequency and is associated with the dissipative effects. The transient mode is responsible for disruption of hydrodynamization process and is controlled with the relaxation time parameter which sets the duration for which viscous effects remain active. These are called the quasi-normal modes of out-of-equilibrium hydrodynamics, analogous to the normal modes of oscillatory systems in classical mechanics.
Right after the collision of heavy-ions, we have a nonequilibrium system of partons for upto 1 fm/c. The fact that applying low-order hydrodynamics does not require local thermalization or even pressure isotropy to show agreement with the measurements [24], had been puzzling, until we discovered that this evolution leads to an attractor [25][26][27][28]. This attractor guides the system evolution to a late time universal trajectory even if initiated with a varied set of starting conditions [29].
The framework of hydrodynamics with initial conditions, followed by a hadron after-burner, has been quite successfully used to explain experimental data obtained from a wide range of systems [30,31]. From most central to ultra peripheral collisions, the system size decreases monotonically. For a constant collisional energy, there should to be a system size below which the QGP droplet will cease to hydrodynamize [32]. Aleksi Kurkela et al., [33,34] has performed a flow analysis with kinetic theory leading to hydrodynamization, through a dimensionless physical quantity called opacity(γ) -a measure of transverse system size in units of the mean free path. As the opacity varies from 0 to 5, the system goes through 3 stages in this order:(a) non-QGP (particle-like) stage, (b) intermediate transition stage and (c) QGP (hydro-like) stage. Ulrich Heinz and Moreland [35] have emphasized considering the multiplicity rapidity density of charged particle -dN/dy along with HBT radii to quantify the smallest QGP size. According to Romatschke [23], the large p T regime of flow is due to non-hydrodynamic mode and this mode can be studied through the relaxation time approach. He suggested that large deviation of elliptic flow (v 2 ) for a variation in shear relaxation time for lowering multiplicity could potentially indicate breakdown of low-order hydrodynamics. The last two of the above studies came to the conclusion that this limit should be around or below dN ch /dy ≈ 2.
The role of relaxation time has been previously analyzed for different settings in hydrodynamics studies [36- 41], including spatial and momentum eccentricity, entropy and elliptic flow for varying relaxation times. However the primary focus of these studies was to find the range of τ π and other second order transport coefficients for which the observables were insensitive, which in turn meant that the magnitude of the second order gradient terms are smaller in comparison to those of first order gradient. In the present work, we check the sensitivity of observables to shear relaxation time in ultra peripheral collision systems to test the breakdown of low-order hydrodynamics. In Sec. II we discuss the framework of the model used. In Sec. II A, we state the initial condition and input parameters involved in the model. Sec. II B describes the observables obtained along with the experimental results in order to fix the centrality related parameters. In Sec. III, we present results of elliptic flow as a function of transverse momentum and multiplicity rapidity density. And in Sec. IV, inferences are drawn based on results obtained along with the possible improvement to this work.

II. FORMALISM
Hydrodynamics is the collective dynamical evolution of a suitably sized bulk medium adhering to the system's symmetries. For the relativistic case, the conservation laws take the form, ∂ µ T µν = 0 for energy-momentum tensor and ∂ µ N µ = 0 for conserved charge. The local values of temperature, T (x), fluid velocity, u µ (x) and chemi- cal potential, µ(x) are chosen as hydrodynamic variables. For ultra-relativistic collisions, where a negligible amount of participating nucleons survive, the conservation equation for baryon number (∂ µ N µ = 0) can be ignored. The energy-momentum tensor can be decomposed as [46]; Here, (energy density) and P (pressure) are scalar coefficients. w µ represents the transverse vector coefficient. ∆ µν ≡ g µν +u µ u ν is the projector operator orthogonal to the fluid velocity(u µ ) and g µν is the space-time metric. The above expression without the Π µν term corresponds to 0 th order ideal hydrodynamics. The Π µν tensor is introduced to account for the dissipative effects and is further decomposed as: Π and π µν are the bulk and shear part of the viscous stress tensor. The form of the shear stress tensor(π µν ) and bulk pressure(Π) are set up in accordance with the covariant form of the second law of thermodynamics [17]. When we set entropy 4−current expression as s µ = su µ , where s is entropy density, we get; where η(shear viscosity) and ζ(bulk viscosity) are the transport coefficients. σ µν (shear tensor) is a traceless, transverse and symmetric tensor. This form of π µν and Π leads to the 1 st order, Navier-Stokes theory. When we introduce perturbations in energy density and fluid velocity, and evolve them, the diffusion speed obtained from the dispersion relation has a form that can increase arbitrarily. This theoretical formulation cannot be considered as a satisfactory one, if it violates causality. It turns out that if the term (−τ π u α ∂ α π µν ) is added in the expression of π µν above, the resulting diffusion speed stays below the speed of light. The coefficient of this newly added term, τ π is called relaxation time. But this is still a makeshift way to restore causality in the system. A good 2 nd order viscous hydrodynamics theory at the very least should reduce to the Navier-Stokes equation in the limit of long wavelengths, and must show causal signal propagation.
Müller [47], Israel and Stewart [48,49](MIS) suggested modification of the entropy 4−current expression used above to include the following term with a viscous stress tensor: where β 0 and β 2 are scalar coefficients. When we use this entropy 4−current in covariant 2 nd law of thermodynamics, the dissipative terms of energy momentum tensor take the following forms [17]: FIG. 3. Charged particle multiplicity rapidity spectra generated (lines) for Au-Au 200 GeV (above) and Pb-Pb 2.76 TeV (below) as a function of number of participants compared with corresponding PHENIX [42] and ALICE [44,45] experimental data (errorbars). The generated data points are labelled with the midpoint of the centrality range in blue color.
Where, ∇ µ = ∆ αµ ∂ α and ∇ α u β is a symbol to represent traceless symmetrization of ∇ α u β . A perturbative analysis with these newly obtained expressions leads to an inherently causal system. There are a few variants of this theory [50], depending on how many terms are kept in π µν and Π expression. The viscous hydrodynamics code used for this study is based on MIS theory. BRSSS theory [51] is a more comprehensive version of MIS hydrodynamics. A few 3 rd order versions have also been worked up [52,53].
The second order viscous hydrodynamics used for this study is a publicly available code 2 , ECHO-QGP [57,58], based on MIS theory. It could be used in either (2 + 1)-D or (3 + 1)-D settings and has been utilized for bulk medium evolution in quarkonia suppression study [59]. Spacetime evolution of all T µν components could be extracted at the output.
2 http://theory.fi.infn.it/echoqgp/index.php A tabular lattice QCD equation of state by Wuppertal-Budapest collaboration [60] has been utilized. In this equation of state, the values for energy density( ), speed of sound(c s ) and pressure(P ) are available starting with the temperature of 100 MeV. In order to get values below this temperature we spline interpolated temperature dependencies of quantities mentioned above with the corresponding values from hadron resonance gas model [61]. Dissipative corrections to the energy momentum tensor in ECHO-QGP are introduced in the same way as stated in Eq. (5). Here the evolution of shear part of the viscous stress tensor is given by [57]; Here, λ 0 is a scalar coefficient and Ω is a traceless, antisymmetric, transverse vorticity tensor. d µ is the covariant derivative given by are the Christoffel symbols. D = u µ d µ , is the comoving time derivative. The evolution of the bulk part of viscous stress tensor is given by; The values of the transport coefficients, τ Π , λ 0 , τ π , η, ζ are required for solving the above two equations, which is obtained from microscopic theory approach to hydrodynamics. τ Π is the bulk viscosity relaxation time, which represents how quickly the above 2 nd order form of bulk pressure relaxes to its leading-order form in Eq.(6). The above two equations are derived under the metric signature choice of (−1, +1, +1, +1).

A. Input parameters
The form of relaxation time has been worked out for hydrodynamics beginning from numerous microscopic theories e.g., Boltzmann theory in the relativistic limit [49,62], weakly coupled QCD [63] and AdS/CFT [51,64,65]. In ECHO-QGP, the relaxation time is introduced as; The coefficient, τ coe here controls the magnitude of shear relaxation time in viscous hydrodynamics. In Sec. (III), we see the consequence of varying this parameter on elliptic flow coefficients for Pb-Pb and Au-Au collisions. A transverse distribution of participating nucleons could serve as an initial condition for hydrodynamics. ECHO-QGP has an optical Glauber model as its default initial condition which assumes independent linear trajectories of nucleons in nuclei that are distributed according to Wood-Saxon distribution [66,67]. Wood-Saxon distribution has a smooth plateau for the nucleus which decays softly towards the edges. Even though the Glauber model does not involve early stage dynamics and fluctuations of any kind, it is still a good approximation nonetheless.
IPGlasma [68,69] is a more realistic initial condition that includes the dynamics beginning from the moment of collision. It is based on the color glass condensate framework. The wavefunction of a nucleus or hadron at high energy could be explained with the effective theory of color glass condensate [70,71]. In IPGlasma model, the color charges inside the nucleons are Gaussian sampled and are taken as the source for gluon fields, which are then evolved using classical Yang-Mills equations [68].  [23]. Phenomenological studies that make use of viscous hydrodynamics have been able to explain flow experimental data only in low pT range. Beyond pT ≈ 4 GeV, presence of non-hydrodynamic mode has been suggested.
We have used the publicly available 3 IPGlasma model that describes a boost invariant (2+1)-D initial state. The energy density in the transverse plane at τ start = 0.2 fm/c for Pb-Pb collision and τ start = 0.6 fm/c for Au-Au collision has been taken as an input for ECHO-QGP. Fig.(1) shows the initial energy densities for 14 centralities as a function of transverse coordinates for Au-Au collision. We ran the [IPGlasma initial condition + ECHO-QGP hydrodynamics] framework for 14 centrality values, with more values near peripheral collisions. The distribution of the nucleons in the nucleus and the distribution of color charge inside nucleons are the FIG. 7. Pion(π + ) elliptic flow coefficient(v2) as a function of transverse momentum(pT ) for 14 centrality classes for Au-Au 200 GeV system obtained with (IPGlasma+2Dhydro) set up along with experimentally measured elliptic flw (blue) from PHENIX [55] for the relaxation/non-hydrodynamic mode decay time, τπ = 3η/sT (green) and 12η/sT (red). The shaded area (yellow) highlights the difference in flow due to variation in relaxation time.
key sources of initial state fluctuations in each collision event. Observables in collider experiments are averaged over a large number of collision events, to account for this event-by-event fluctuation. For both Au-Au and Pb-Pb collision systems, we produce an initial state with 400 different sets of nucleon positions that are then combined into one. The total inelastic nucleon-nucleon cross section is set to 61.8 mb for Pb-Pb system and 42 mb for Au-Au system in both, IPGlasma and in hydrodynamics, taken from Monte Carlo Glauber analysis [72]. Shear viscosity to entropy density ratio(η/s) is taken as a constant, 0.1 (≈ 1.25× 1 4π ) [3], which is above the theoretical minimum KSS limit [73]. Bulk viscosity has not been included in this study. The pseudo-critical temperature, at which quarks to hadron phase transition occur, has been calculated by various lattice QCD collaborations, is an input parameter. It is set to the recently calculated value of 156 MeV [74]. Chemical freezeout is a point at which the inelastic scatterings cease to exist between produced hadrons. This point is decided by the temperature, which in the present model is fixed at 150 MeV [75]. Fig. 2 shows the p T spectra of pions(π + ) produced for the two mentioned collision systems along with corresponding experimentally measured p T spectra. The generated spectra adequately comply with experimental values only in low p T regime, where the hydrodynamic mode operates. The energy density profile plotted as a function of transverse coordinate from IPGlasma had to be scaled before being used in hydrodynamics. Fig. 1 shows this scaled energy density distribution. This fixed the energy density scaling parameter such that the produced p T spectra and the maxima of rapidity spectra(dN/dy) at each centrality matches with the corresponding experimental measured data for both collision systems. Fig. 3 shows rapidity spectra normalized to N part /2 as a function of N part . In addition to energy density scaling, the rapidity spectra had to be scaled to match with experimental results as shown in Fig. 3. For the Au-Au system, charged particle normalized rapidity spectra were scaled up by a factor of 2, whereas for Pb-Pb system, this scaling was 6, and the corresponding scaling used for pions was 3.6. We chose a centrality range spaced by 5% in peripheral collisions except for the last centrality class, 90-100%. The impact parameter and N part values for all of these centrality ranges are taken from a Monte Carlo Glauber analysis [72]. The reason for taking more values towards the peripheral side was to capture fine variations of flow for decreasing dN/dy as could be seen in Fig. 10 in Sec. III. However there was no experimental reference to set parameters for these in-between centrality values for p T -spectra and dN/dy vs N part plot. Hence, we selected two values around each experimental centrality point starting from 60% as could be seen in dN/dy vs N part plot (Fig. 3). There was no experimental point at 90-100% so we settled with just one extrapolated value which follows the trend of data. The blue labels on data points in Fig. 3 are the mid centrality value of that data point. For calculating observables for charged particles, we have added the corresponding values for the pions(π + +π − ), kaons(K + +K − ) and protons(p + +p − ) since these are abundantly produced species in high energy collisions. Momentum space eccentricity which is the precursor of elliptic flow can be calculated in terms of T µν components as:   FIG. 9. pT integrated elliptic flow in proton-proton collision at 7 TeV produced using SONIC model, as a function of multiplicity pseudorapidity spectra for the mentioned values of η/s and ζ/s. For η/s = 0.08 and ζ/s = 0 (blue), elliptic flow has errorbar due to variation in shear relaxation or nonhydro mode decay time, which increases in size for decreasing dN/dη. Plot taken from [23].

B. Fixing centrality parameters
ECHO-QGP calculates this quantity for ideal hydrodynamic case, which takes the form: To generate momentum eccentricity for the viscous case, we have modified the above expression by adding viscous component term, (π xx + π yy ) to the integrand in both numerator and denominator. Fig. 4 shows spatial eccentricity( c ) and momentum space eccentricity( p ) for Au-Au and Pb-Pb collision, generated at 50 − 60% centrality for the two mentioned shear relaxation times. Momentum anisotropy quantified by momentum eccentricity increases at the expense of spatial anisotropy quantified by spatial eccentricities along the evolution [16]. The variation in nonhydrodynamic mode decay time seems to have negligible effect on spatial eccentricity. The distinguishing feature between the two systems is that the early time c for Pb-Pb decreases more rapidly than that for Au-Au collisions. Below the pseudo-critical temperature, hadronic picture should emerge. Particles of various species are assigned momentum according to Cooper-Frye scheme [76]. The resulting momentum spectrum is then used to calculate the elliptic flow, v 2 = cos[2(φ − Ψ RP )] , where Ψ RP is the reaction plane angle which acts as a reference plane and φ is the transverse plane angle for a given particle with respect to the reaction plane. Fig. 5 shows the average transverse momentum evolution as a function of centrality. Results for the two values of shear relaxation time have been plotted and compared with experimental values for pions. We notice, that the model show agreement with experimental values for most of the centrality classes apart from the peripheral ones. The values for Pb-Pb collisions had to be scaled up by a factor of 1.3. This could be due to underproduction of hadrons in the hydrodynamics, because the multiplicity has been used as the weight factor for calculating mean p T .

III. FLOW RESULTS AND DISCUSSION
Romatschke [23] has put forth a quantitative test for applicability of hydrodynamics by checking the sensitivity of certain observables(like elliptic flow) to the nonhydrodynamic mode. The idea is that hydrodynamics can be used to describe a system if the non-hydrodynamic mode is sub-dominant and there exists a local rest frame. With QCD as the microscopic theory, approximate transverse momentum range of hydrodynamic mode is 3 to 7 GeV. Fig. 6 illustrates this p T -range where hydro and non-hydro modes operate. This is what we have tried checking for Au-Au 200 GeV in Fig. 7 and for Pb-Pb 2.76 TeV in Fig. 8. Peripheral collisions are the system of interest, but experimentally measured anisotropic flow results are only available upto 50-60% centrality class. We hence presented the results for the complete centrality range. In Fig. 7, for 0-5%, 5-10% and 10-20% centralities, we see no separation between elliptic flow curves for non-hydrodynamic mode decay times, τ π = 3η/sT and 12η/sT . From 20-30% centrality class onwards we notice the separation between these two flow curves to be increasing. Experimental data has been plotted just for reference that show our results are quite close to experimentally measured flow results. The important point to notice is that along increasing centrality, the point at which the two flow curves separate shift towards lower p T values. Which means that with increasing centrality and decreasing system size, the hydrodynamic mode is shrinking and non-hydrodynamic mode is getting dominant. Hence in a way we are witnessing limit of applicability of low-order hydrodynamics for decreasing system size at constant collisional energy(here, 200 GeV). Fig. 8 shows p T dependence of pion(π + ) elliptic flow with complete centrality range for τ π = 3η/sT and 8η/sT . We notice all the structures mentioned above for Au-Au, 200 GeV system. We notice a better match between produced elliptic flow and experimental data 10-20% onwards. We chose pions for this analysis because they are the lightest of particle species produced and hence adequately represents the bulk medium. One additional point to notice is, for 10 − 20% centrality in Au-Au collisions and classes 0−5%, 5−10% in Pb-Pb collision system, our model fails to reproduce the measured This abrupt change in elliptic flow is indicative of breakdown of hydrodynamics, and it is seemingly happening at roughly dN/dη < 2 in Fig. 9. We tried checking this feature in our (IPGlasma+2Dhydro) analysis as shown in Figs. 11 and 10. Fig. 10 presents the un-normalized p T integrated elliptic flow as a function of multiplicity rapidity density(dN/dy). The data points from our analysis are labelled by the centrality class in order to track the point at which flow changes abruptly between the relaxation time curves. This is why we selected more centrality points in peripheral collision side. We notice a steady increase in separation between the two flow curves for both Au-Au and Pb-Pb system which is in reasonably close agreement with Romatschke's work.
We also notice that the two relaxation time flow curves of same centrality do not have same multiplicity rapidity density value (the x co-ordinate). This would mean, that for an increase in relaxation time, flow shifts to a lower multiplicity value. We also notice that the flow for τ π = 3η/sT for both, Au-Au and Pb-Pb systems, acquire negative values, which is also apparent from the elliptic flow for 90-100% centrality class in Figs. 8 and 7. Fig. 11 shows normalized p T integrated elliptic flow as a function of charged particle multiplicity rapidity density(dN/dy) for peripheral collisions. We clearly notice the sudden increase in separation of flow curves for the two mentioned relaxation times for both the collision systems. But we don't have a centrality resolution good enough to decide the onset of hydrodynamization. A approximate limit we can deduce from Fig. 11 is dN/dy ≈ 10 which is quite larger than the prediction of dN/dη < 2 [24,35]. However if the hadron resonance gas to de-confined quarks transition in high temperature regime is a crossover, we expect to find a region where analysis would be indecisive like what Aleksi Kurkela et al. obtained [33,34]. The problem lies in the absence of experimental reference data to set the scaling parameter of IPGlasma for such high centrality classes.

IV. CONCLUSION AND OUTLOOK
In this study we analyze the non-hydrodynamic mode in an attempt to find the onset of hydrodynamization in peripheral collision system of Au-Au and Pb-Pb at 200 GeV and 2.76 TeV center of mass per energy nucleon, respectively. We use the energy density profile from color glass condensate based IPGlasma model as the initial condition in 2D ECHO-QGP which is a 2 nd order viscous hydrodynamic code based on MIS theory. p T spectra and multiplicity rapidity density(dN/dy/(N part /2)) as a function of N part is used to constrain the centrality scaling parameter of IPGlasma. Mean p T as a function of centrality, evolution of spatial and momentum eccentricity has also been generated for both the systems. The shear viscosity to entropy density ratio is set as η/s = 0.1 and bulk viscosity has not been considered in this work. We study the variation in the strength of nonhydrodynamic mode through the shear relaxation time, whose value is set to (3 − 12)η/sT for Au-Au system and (3 − 8)η/sT for Pb-Pb system. Elliptic flow generated as a function of p T is compared with 2 nd anisotropic flow coefficient from experiments for the above respective values of relaxation time, for all of the 14 centrality classes. Normalized and un-normalized p T integrated elliptic flow has been studied as a function of multiplicity rapidity density in peripheral collisions especially. We found the following: • From p T dependence of elliptic flow across centralities for Au-Au in Fig. 7 and for Pb-Pb in Fig.  8, we found that the shear relaxation time does control the non-hydrodynamic mode of the system as predicted by P. Romatschke. This inference was guided by the observation that the point after which the flow for the two relaxation times separate sharply from each other, shifts to lower p T values for increasing centrality classes (or decreasing system size at a constant energy of collision).
• We later attempted testing the onset of hydrodynamization from charged particle multiplicity rapidity density dependence of p T integrated elliptic flow. We did notice an abrupt increase in flow for decreasing system size or number of par-ticipants, indicating increased dominance of nonhydrodynamic mode and simultaneous breakdown of hydrodynamic description. However we could not resolve the dN/dy below the value of 10 enough to quantitatively decide the onset point.
• We found a good agreement between the generated p T dependence of elliptic flow results and the measured flow data from PHENIX and ALICE Collaborations for Au-Au and Pb-Pb systems, respectively, except near most centrality of 10-20% class for Au-Au collisions and and of 0-5% and 5-10% class for Pb-Pb collisions.
There is significant scope for improving this framework further by including an after-burner stage that will incorporate hadron resonance decays and scattering which could affect the generated flow [77]. It will be interesting to compare lowest fluid size from other methods in the future work. The initial state involvement could also be improved by using more components of T µν in hydrodynamics [78,79]. One can also switch to 3−D IPGlasma initial condition [80]. Bulk viscosity has been kept zero in this study. But it does play significant role in evolution [81]. The relaxation times for bulk viscosity could be independently analyzed. It will be interesting to see if elliptic flows of different particle species diverge for decreasing rapidity spectra at different points. If they do so, it would support the idea of multiple-fluid scenario in heavy ion collision. This study could be extended to small and lower energy system where the netbaryon potential is non-zero, for which particle current conservation should be included [82]. In a recent study, Plumberg et al. [83], have conducted a causality analysis of each fluid cell of hydrodynamics for its complete evolution. They found causality being violated of nonhyperbolic(v 2 < 0) and superluminal(v 2 > c 2 ) type at early times in evolution. This violation is significantly reduced if a pre-equilibrium stage like KøMPøST [84] is used. It will be interesting to see the repercussions of such a study on onset of hydrodynamization.

ACKNOWLEDGMENTS
We are thankful to Gabriele Inghirami for clearing our doubts and helping at numerous times in using the viscous hydrodynamics code. We are also grateful to Chun Shen, Paul Romatschke and Rajeev Bhalerao for clearing our doubts about flow. Nikhil Hatwar would like to thank BITS-Pilani for the financial support.